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Summary. A mayor problem that arises in the computation of stellar atmosphere 
models is the self consistent determination of the temperature distribution via the 
constraint of energy conservation. The energy balance includes the gains due to the 
absorption of radiation: J xii^) J {'^)di', SLiid the losses due to emission: J x{i^)S{u)du . 
It is well known that within each one of the two above integrals the part correspond- 
ing to spectral ranges whose opacity xi'^) is huge can overcome by many orders of 
magnitude the part that corresponds to the remaining frequencies. On the other 
hand, at those frequencies where xi'^) is very large, the mean intensity J{i^) of the 
radiation field shall be equal, up to many significant digits, to the source function 
S{h') and consequently to the Planck function B{v,T). Then their net share to the 
energy balance shall be null, albeit separately their contributions to the gain and 
loss integrals are the most important numerically. Thus the spectral range whose 
physical contribution to the overall balance is null will dominate numerically both 
sides of the energy balance equation, and consequently the errors on the determi- 
nation of Jiy) and Siy) at these frequencies will falsify the balance. It is possible 
to circumvent the numerical problem brought about by the foregoing circumstances 
by solving the radiative transfer equation for the variable /(n, v) — S{v), instead of 
the customary intensity /(n, v). 

We present here a novel iterative algorithm, based on iteration factors already 
employed by us with success, which makes it possible a fast correction of the tem- 
perature by computing directly the difference between the radiative losses and gains 
at each step of the iterations. 



1 A Two-Block Iterative Scheme for Modeling Stellar 
Atmospheres 

The distinctive feature of our algorithm for an iterative temperature correction 
procedure comes from the organization into two fundamental macro-blocks of 
the ensemble of the mathematical equations which represent the physical phe- 
nomena that occur in a stellar atmosphere. Within each block a self-consistent 
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physical problem is solved by assuming that the current values of some input 
variables arc known, so that the updated values of the output variables of the 
block can be obtained. 

The first macro-block (the structural one) consists of the equations that 
express the conservation of mass and momentum. The latter, together with 
the equation of state, are the backbone of the "hard " structure of the model. 
If we assume that the temperature distribution T(r) is given, we can deal 
with the equations of this block in a quite easy way, and get the distribution 
of the pressure P(r) and density p(r), hence the values of all the relevant 
thermodynamical and transport coefficients. All these values depend of course 
on the initial temperature distribution assumed. 

Once we know the hard structure, namely P{r) and p{r), as well as the 
ensemble of the thermodynamical and transport coefficients at each point of 
the physical system, we can tackle the radiative transfer (RT) equations, cou- 
pled through the constraint of the conservation of energy, in order to obtain 
with relative ease a new temperature distribution. These equations, together 
with the energy conservation constraint, form the second macro-block (the 
energetic one). The two macro- blocks are coupled both physically and mathe- 
matically via the equation of state. However this coupling is not very strong, so 
each block can be treated mathematically in an independent way by means of 
a sequential procedure that is iterated till the required numerical convergence. 
According to our experience, such an iterative procedure quickly converges (in 
about ten iterations) , and numerical difficulties never show up. An exhaustive 
presentation of the algorithm can be found in [1]. 

2 The Iterative Scheme for the Temperature Correction 

The aim of this work is to discuss the numerical treatment of the second 
macro-block, given the hard structure of the system. Such a procedure is usu- 
ally called the temperature correction, because it is matter of determining - 
also iteratively - the temperature distribution that satisfies simultaneously 
both the RT equations and the relevant constraint of total energy conserva- 
tion, through the correction of a previous temperature distribution T *(r), 
assumed to be close to the required one. For the sake of an easier presenta- 
tion, we will consider here the instance, in which energy is transported only 
by radiation. It is simple to generalize the algorithm in case that both ra- 
diative and convective transport are important, as shown in [2]. Likewise, we 
will consider the paradigm case of a model atmosphere under the simplifying 
hypotheses of hydrostatic and local thermodynamical equilibrium. 

The RT process is described by means of a set of kin(^tic equations, one 
for each of the specific intensities /(n, i/) employed for the representation 
of the radiation field. The corresponding sink terms, i.e. the specific radia- 
tive energy removed at each point, are assumed to be proportional to /(n, v) 
through the absorption coefHcient a(i/) and the diffusion coefficient <j(i/) . The 
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source terms are described by introducing the thermal emission r]{i/), which 
depends at each point on the local temperature according to the Kirchhoff's 
law: r}{i') = a{i')B{i', T), and a diffusion (scattering) term a{i')J{i'). The sym- 
bol J(i/) denotes the monochromatic mean intensity of the radiation field, i.e. 
the integral of /(n, i/) over all the directions n times l/Air. In the processes of 
scattering the amount of monochromatic radiative energy gained compensates 
exactly the energy lost, when all the directions are taken into account. 

Under the assumption of radiative equilibrium (RE), the constraint of 
energy balance reads 



which expresses the equality between thermal emission and the absorption 
processes. The rhs of eq. (1) is straightforwardly computed from the solution 
of the RT equations. If this term is known, eq. (1) becomes an implicit tran- 
scendental equation that yields at each point of the atmosphere the local value 
of T{r) in a trivial way. It must be stressed that at this stage of the procedure 
the absorption coefficient a(z^) is assumed to be known. The only difficulty is 
that J{v) comes from the solution of the RT equations, which depends on the 
temperature through the Planck function B{i', T) that enters into the relative 
source terms. 

This dependence entails introducing a new iterative procedure, which 
can be sketched as follows. We start from a trial temperature distribution 
T * = T *(r), which allows us to solve explicitly the RT equations, namely to 
get the specific intensities J(n, v) at each point. We can then compute the cor- 
responding mean intensity J ^(v), from which we obtain the explicit value of 
the rhs of eq. (1), corresponding to the current value T * of the temperature. 
By solving eq. (1) as a transcendental equation we get the new tempera- 
ture distribution T Unfortunately the rate of convergence of this simple 
scheme, usually called the yl-iteration, is exceedingly slow; in the numerical 
practice often it does not converge at all. 

Although the yl-iteration does not converge, the relevant working plan is 
so simple that it deserves to be kept. A slight alteration of the above itera- 
tive procedure can dramatically improve the rate of convergence. Moreover - 
which is even more remarkable - that neither implies a heavy computational 
burden nor introduces the risk of numerical instabilities. The change has to 
be introduced once J ^{v), the frequency integrated mean intensity J ' and 
the total energy absorbed J^^ have been computed from the RT equations. 
The improvement consists in the re-computation of the total energy absorbed, 
whose value at each point we will denote now by °. These values are eval- 
uated straightforwardly from the system of moments of the RT eqiiation, and 
computed in such a way that they meet certain properties of the converged 
solution. These properties, relevant to radiative transfer, are a consequence of 
the constancy of the total flux H, i.e. the integral over all the frequencies of 
the monochromatic flux H(y). 




(1) 
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Let us consider in detail this new intermediate step. The RT equation for 
a one-dimensional planar medium is 

11 - x(j^) [I{^^,>^) - s{i^)] , (2) 

where ji is the cosine of the angle between the direction of propagation n 
and the r-axis, and x(z/) = a{v) + a{v) the total opacity. The source function 
is expressed as 

S{u) = e{v)B(v,T) + [1 - e{v)] J{v) , (3) 
where e{u) is defined as 

The equations for the three first y^-moments of the specific intensity: 
J(z^), H{i') and K{v), are obtained by integrating eq. (2) over d/x and /xd/x 
respectively: 



^ _a{v) [J{v) - B{v,T)] (5) 



dr 



and 



dKiv) , •, s ,„s 



dr 

By integration over all the frequencies, we get the bolometric equations 



^ = B,- J, (7) and ^ = - , (8) 

dr dr 



where J, H and K are the frequency-integrated zero- , first- and second- 
order /i- moments of the specific intensity. The quantities Ba, Ja and H^. are 
the integrals over all the frequencies of B{v,T) and J{v) weighted with the 
absorption coefficient a{v)^ and that of H{y) weighted with the opacity x(^^)) 
respectively. 

We recognize in eq. (7) the statement of the energy balance: the difference 
between the total emission and the total absorption is compensated by the 
variation of the total radiative flux. Under the condition of RE, this difference 
must be zero, therefore the variation of H will be null. The constant total 
radiative flux H* that corresponds to the situation of RE is directly related 
to the luminosity, viz the effective temperature Tg// of the star: 



H* = ^ Ttff , (9) 



47r 

where ur = 5.6696 • 10~^ erg ■ cm~'^ ■ ■ is the Stefan-Boltzmann 
constant. 
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Once introduced these preliminary definitions, we can get back to the 
iterative correction of the temperature. In the first part of each step of iter- 
ations, given the current value of the temperature T * and consequently that 
of B{i>,T '), we can compute not only J '(z/), but also H and K 
hence, trivially, the frequency- integrated quantities: J % J^, % and 
K The extra computational work for these quadratures adds very little to 
that required for the quadrature of J{v). Then, by making use of the foregoing 
quantities, we compute at each point the following ratios: 

f (10) , /? = (11) and aj ^^f. (12) . 

Because these ratios are the quotient between homologous quantities (in 
the practice pairs of quantities with the same physical behavior), they shall 
result almost independent of the input B{v,T *). We can consider them as 
quasi-invariant, and employ them as iteration factors. 

Afterward we have obtained - in such an easy way - the iteration factors, we 
will proceed recomputing the improved value of the energy absorbed: 
Thus, by recalling the definition of the factor /? given by eq. (11), we can 
rewrite eq. (8) in the form 

f = - ■ 

where we have introduced the value of H* , a data that characterizes the 
problem under consideration, as shown by eq. (9). The quasi-invariant itera- 
tion factor /? is almost independent of the input B{y,T *). That allows us to 
obtain an improved value K ' for the moment which takes into account 
the correct value of the radiative flux. The first iteration factor /, defined by 
eq. (10), is the Variable Eddington Factor. It shall yield the corresponding 
improved value of J ' 

The third quasi-invariant factor a,/, defined by eq. (12), shall furnish us 
the up-dated value of the total energy absorbed, i.e. = aj J ^ This 

is the rhs of the equation that expresses the constraint of RE. We will solve 
the transcendental equation 

a{v)B{v,T '+!) dv = " , (14) 

by using now the improved value J^' °, in order to get the up-dated tem- 
perature T at each point. 

On one hand the iteration factors are almost independent of the input 
B{i>,T *); on the other they are very easily computed from the RT solution. 
Thus the computation of the improved value J^^ of the total energy absorbed 
is as simple in the practice as the previous computation of J^^ directly from the 
RT solution. But now J,j* includes the information on the constancy of the 
total radiative fiux H*, namely the condition of RE: by employing the iteration 
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factors we solve the set of the RT equations consistently with the constraint of 
energy conservation. That makes the iterative procedure extraordinarily fast: 
the rate of convergence improves dramatically at the price of a little extra 
computational cost. 

3 A Pending Difficulty with the Temperature Correction 

It is well known , however, that within each one of the integrals in eq. (1) the 
contributions due to the spectral ranges corresponding to the most opaque ra- 
diative transitions can overwhelm by many orders of magnitude those due to 
the remaining frequencies. On the other hand, far from the superficial layers 
the monochromatic mean intensity J{i^) will be equal up to many significant 
digits to the source function S{h'), and consequently to the Planck function 
B{v,T), at those frequencies for which the absorption coefficient a(y) - hence 
the opacity x(i^) - is very large. Then the net contribution of the latter to the 
energy balance must be null, albeit separately their contribution to the two 
integrals are far the most important numerically. Thus those spectral ranges, 
whose contribution to the overall balance is null, will dominate numerically 
both side of the relevant equation, and consequently the errors on the deter- 
mination of J(i/) and S{y) at these frequencies shall falsify the balance. 

A way to circumvent this severe numerical problem is to solve the RT 
equations in the new variable /(n, v) — S{v)^ instead of the specific intensity 
/(n, v) customarily employed. Then it will be possible to compute the dif- 
ferences J{u) — B{v,T) directly from the solution of the RT equations. The 
spurious contributions above mentioned disappear, because the differences au- 
tomatically vanish when it is the case, independently of the absolute vales of 
/(n, u) and B{v, T). The introduction of the new variable /(n, i') — S{v) comes 
in a natural way from the use of our Implicit Integral Method (see [3] and [4]). 
Its application to the problem under study here have been preliminary dis- 
cussed in [5], both from the physical and the mathematical standpoint. 

By combining the above way of solving the RT equations with the method 
of the iteration factors, it is possible to correct the temperature in spite of the 
difficulties mentioned previously. In the first part of each step of iterations, 
given the input B{v,T '^), we will compute directly the differences J{v) — 
B{v, T *) for each frequency. Hence, by adding to the latter the data B{u, T 
we will obtain the monochromatic mean intensities J{v), as well as all the 
other moments, both monochromatic and frequency-integrated, that we need 
in order to compute the iteration factors defined by eqs. (10) through (12). 
It must be stressed that we get directly the differences J{v) — B{v,T *), and 
consequently their product by the absorption coeflacient a{v) integrated over 
all the frequencies, i.e. 

(J - B): = / a{u) [Jiu) - B{u,T')] du , (15) 
Jo 
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is evaluated without the need of computing separately its two components 

and B^^. If (J — B)^ is null at each depth point, the condition of RE is 
fulfilled everywhere, and therefore the process of correcting the temperature 
is over. If not, we must proceed with the second part of the iteration step to 
determine a new temperature distribution T via an equation like eq. (15). 

Once we have obtained the values of (J — B)^ on the one side, and the 
monochromatic and frequency- integrated moments (in particular J J) as well 
as the iteration factors on the other, we go forward to compute the improved 
values J * and J^* °. As already said before, on one hand the dependence of 
the latter on the input B{u, T ') is weaker than that of J^^, on the other they 
take into account the information on the constancy of the total radiative flux 
in a stronger way. The effect brought about by these advantages is measured 
by the factor 

Tie Tic 

whose value converges very quickly to unit. 

Let us get back now to cq. (15), which is the protagonist of the second part 
of each step of iterations. We can add to each of his members the quantity 
"fBaiT '), where Ba{T *) is the input for the current step of iterations, in 
order to get 

/ a{v) [B{u,T'+') - jB{u,T')] du = Jj - lB^{T = 
Jo 

= 7 [Jj - B^T')] = 7(J - B): . (17) 

The term {J — B)^, defined by cq. (15), is actually what we compute di- 
rectly as a whole, starting from the solution /(n, i^) — S{i') of the RT equations. 
In such a way there vanish automatically the differences between absorption 
and emission for those frequencies, whose net contribution to the energy bal- 
ance must be null, in spite of the fact that the corresponding values of a{u)J{v) 
and a{v)B{u, T) are overwhelming (sometimes by many orders of magnitude) 
in the corresponding integrals. 

That was exactly our aim. The solution of the new transcendental equation 
(eq. [17]), that we have to solve now, presents the same degree of difficulty as 
the previous equation (14). We can conclude that to iterate by making use of 
the iteration factors not only speeds up dramatically the rate of convergence 
of the procedure, but also allows us to take into in account in the numerical 
algorithm only the intrinsic effects of the physical process that are actually 
necessary. 

The term (J — B)^ represents a local correction: at those points where it 

is different from zero the condition of RE is not satisfied, and consequently 
it will be necessary to correct the temperature there. However it may occur 
that (J — B)^ be null at a certain point, and nevertheless the correction be 



8 O. Cardona et al. 



necessary, because the condition of RE is not fulfilled elsewhere. This non-local 
effect is taken into account (iteratively) through the factor 7. 

4 Conclusions 

We have presented here a novel method for building stellar atmosphere mod- 
els, more; precisely, for the calculations required by the determination of the 
temperature at each point. Because of radiative transfer one must take into 
account the coupling among the physical conditions at each and every point of 
the system. In particular, the temperature at each point is linked with those 
at all the others. We ideally remove the coupling, and proceed sequentially to 
get the physical conditions at each point under the assumption that those at 
all the other points be known. In order to achieve the consistency - i.e. the 
coupling - we developed an iterative procedure. Such a treatment constitutes 
a numerical simulation of the physical processes that are at the origin of the 
corresponding equations. We tackle and solve those problems that arise when 
the fluxes to be determined are the difference between quantities (densities) 
that are very close, so that the difference between their values is less than the 
numerical accuracy of the computations. 

As a proof of the quality of the method proposed, we claim that the in- 
tegrated radiative flux H computed at each step of iterations, quickly con- 
verges (in a few iterations) to the actual value H* . The inaccuracies on the 
converged numerical values of H * are due to the quadrature of H ^{f) over 
frequencies, inaccuracies that we estimate to be more close to 0.01% than to 
0.1%. As the problem with the numerical quadratures over frequencies is the 
same as above, we impose that the equality between the two terms of the tran- 
scendental equation (14) be numerically accurate to one order of magnitude 
less, namely to 0.001%, in order not to increase the unavoidable indetermina- 
tion already brought about by the numerical treatment of the data. 
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